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PREFACE 
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of  ongoing  work  at  the  U.  S.  Army  Engineer  Waterways  Experiment  Station 
(WES),  in  Civil  Works  Investigation  Study:  "Seismic  Effects  of  Reser- 
voir Loading  and  Fluid  Injection,"  sponsored  by  the  Office,  Chief  of 
Engineers,  U.  S.  Army. 

This  study  is  directed  by  Dr.  E.  L.  Krinitzsky,  Engineering  neol- 
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General  direction  was  by  Mr.  J.  P.  Sale,  Chief,  GL,  and  Dr.  D.  C.  Banks, 
Chief,  EC&RMD. 

COL  J.  L.  Cannon,  CE,  and  Mr.  F.  R.  Brown  were  Director  and 
Technical  Director,  respectively,  of  WES  during  the  period  of  this  study. 
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aeismographlc  ■Cations  around  a large  raaarvoir.  In  racant  years,  the  prac- 
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tics  of  placing  such  networks  around  large  reservoirs  to  halp  In  tha  mitiga- 
tion of  sslsmlc  risk  haa  become  widely  accepted.  Tha  purpose  of  the  sensitive 
seismographs  la  a)  to  determine  the  frequency  of  the  local  earthquakes  (if 
any),  b)  to  determine  the  location  of  the  aelsmlc  activity  and  lta  depth, 
c)  to  determine  the  magnitude  and  some  Indication  of  the  focal  mechanisms 
of  tha  earthquakes,  and  d)  to  allow  tha  prediction  of  the  fluctuations  in 
time  of  the  earthquake  occurrence.  It  haa  bean  also  widely  recognised  that 
in  order  to  obtain  a comparative  background.  It  la  essential  to  install  the 
aenaltlve  salamographlc  network  in  the  vicinity  of  the  dam  before  major  con- 
atructlon  begins. 

Although  thia  practice  of  making  precise  aelamlclty  studies  around 
large  reservoirs  haa  become  commonplace.  It  la  remarkable  how  little  atten- 
tion haa  been  given  to  the  optimum  designs  of  the  networks . The  problem  Is 
a particular  case  of  a more  general  one  which  relates  to  tha  modern  practice 
of  hypocentral  location  in  local  regions.  Considerable  progress  has  been 
made  in  this  field  through  the  use  of  algorithms  that  depend  on  high  speed 
computers  (Fllnn,  1960;  Nordqulst,  1964;  Turcotte,  1964). 

The  accuracy  of  the  current  hypocentral  location  algorithms  now  used 
depends  upon  the  adequacy  of  the  seismic  travel  times  (or  velocity  variations 
In  tha  crust)  that  are  adopted.  The  usual  situation  is  that  the  adopted 
values  have  been  determined  prior  to  the  hypocentral  location  attempt,  often 
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from  a set  of  travel-time  observations  made  at  a different  location.  Thus, 
for  the  site  In  question,  such  as  a new  dam,  the  travel-time  model  used  will 
contain  errors  that  may  result  In  grossly  biased  models,  which  in  turn 
produce  algnlf icantly  biased  locations.  It  has  now  been  shown  (Dewey,  1971; 
Douglas  and  Lilwall,  1972;  Douglas,  1967;  Bolt  and  Freedman,  1968)  that  there 
la  a marked  advantage  In  adjusting  the  preliminary  locations  of  earthquakes 
simultaneously  In  groups,  rather  than  one  by  one.  The  method,  called  joint 
hypocenter  determination,  is  a special  application  of  analysis  of  variance  in 
which  arrival-time  measurements  for  many  earthquakes  in  the  region  are  pooled. 
The  result  is  a more  reliable  estimate  of  the  relative  adjustments  between 
Individual  hypocenters  of  the  group.  Systematic  displacements  of  the  group 
as  a whole  and  in  physical  instabilities  in  the  matrix  inversion  are  minimised 
by  the  imposition  of  side  conditions  and  penalty  equations  that  restrict 
adjustments  to  small  steps  at  each  iteration. 

One  proven  procedure  to  obtain  realistic  side  constraints  is  to  include 
in  the  simultaneous  group  reduction  the  observed  times  of  an  earthquake  with 
a prior  well-determined  location.  For  convenience,  the  adopted  parameters 
of  this  master  earthquake  may  be  held  essentially  fixed  by  the  application 
of  appropriate  weights.  The  master  earthquake  may  be,  for  a reservoir,  a 
specially  detonated  chemical  explosion  with  known  origin  time  and  location. 

The  discussion  in  this  report  makes  some  use  of  recent  theoretical 
advances  in  geophysical  inverse  theory  and  indicates  how  in  more  advanced 
research  applications  Inverse  theory  can  treat  the  problem  of  optimum  dis- 
tribution of  selsmographlc  stations  in  a network  and  the  optimum  analysis 
of  arrival  times  of  various  phases  from  earthquakes  from  within  the  network 
(Crosson,  1976).  This  type  of  description  permits  a better  understanding 


of  Che  mathematics  involved  and  provides  measures  of  the  resolution  and 
precision  of  the  calculated  earthquake  parameters  and  crustal  structure. 

The  main  purposes  of  this  report,  however,  are  to  provide  practical  guidelines 
to  engineers  and  geologists  who  wish  to  establish  a well-designed  selsmographic 
network  around  an  engineering  facility  and  to  provide  a tested  and  stable 
computer  program  that  will  allow  straightforward  but  high  quality  earthquake 
locations  to  be  calculated.  For  this  reason,  the  more  mathematical  details 
of  Inverse  theory  in  this  context  are  not  included. 

1.2  Seismic  Risk  to  Dams 

A large  earth  or  concrete  dam  represents  a notably  important  type  of 
seismic  hazard  evaluation  problem.  Not  only  is  the  dam  in  itself  a rela- 
tively expensive  project,  but  it  is  intimately  involved  in  the  whole  economy 
through  power  generation,  flood  control,  irrigation,  etc.  In  addition, 
structural  failure  of  the  dam  may  lead  to  major  disaster,  because  large 
populations  may  be  exposed  to  sudden  flooding.  Major  damage  has  occurred  to 
earth  dams  from  natural  earthquakes;  for  example,  the  Hebgen  Dam  in  Montana 
and  the  Lower  Van  Norman  Dam  in  California.  Many  large  dams  in  the  United 
States  are  located  in  highly  seismic  regions  close  to  areas  that  have  in  the 
past  suffered  from  major  earthquakes.  It  is  therefore  necessary  to  keep  in 
mind  the  likelihood  of  future  damaging  shocks  to  these  structures. 

In  addition  to  the  hazard  from  natural  earthquakes,  several  cases  have 
now  been  documented  in  different  countries  of  damaging  earthquakes  apparently 
related  in  some  way  to  reservoir  loading  behind  the  dam  (Bolt  et  al,  1977). 

Some  of  these  examples  have  occurred  in  regions  that  have  not  been  thought 
to  be  even  moderately  seismic.  The  clearest  cases  are  a)  Lake  Kariba  in 
Central  Africa,  the  world's  largest  artificial  reservoir,  b)  Koyna  in 
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India,  and  c)  Hsingfenkiang  in  che  People's  Republic  of  China.  In  these 
examples,  the  largest  earthquakes  reached  a magnitude  of  6.4.  At  Koyna,  in 
addition  to  significant  cracking  of  the  concrete  gravity  dam  which  required  a 
major  repair  and  strengthening,  numerous  collapsed  houses  in  the  vicinity 
caused  loss  of  life.  Hsingfenkiang  Dam,  which  is  located  in  an  essentially 
aselsmic  region,  impounded  a reservoir  which,  subsequent  to  1959,  was  the 
site  of  numerous  small  shallow  earthquakes.  The  principal  shock  of  the  series 
in  1962  had  a magnitude  6.1  and  produced  a crack  82  meters  long  in  the  upper 
concrete  dam  structure. 

The  cases  of  induced  seismicity  now  make  it  necessary  to  consider  the 
risk  from  Induced  seismicity  for  all  proposed  large  dams. 

The  idea  that  earthquakes  might  be  triggered  by  impounding  surface  water 
is  not  new.  In  the  1870's  the  U.S.  Corps  of  Engineers  rejected  proposals 
for  major  water  storage  in  the  Salton  Sea  in  southern  California  on  the  grounds 
that  such  action  might  cause  earthquakes.  The  first  detailed  evidence  ol  such 
an  effect  came  in  the  United  States  with  the  filling  of  Lake  Mead  behind 
Hoover  Dam  (height  142  meters),  Nevada  , beginning  in  1935.  Although  it  is 
not  certain  that  there  was  no  local  seismicity  before  1935,  it  is  sure  that 
after  1936  earthquakes  were  much  more  common.  A network  of  seismographs  was 
specially  installed  around  Lake  Head  and  showed  that  after  the  largest  earth- 
quake (magnitude  about  5)  in  1940  the  seismicity  has  declined.  Hundreds  of 
detected  earthquakes  cluster  on  steeply  dipping  faults  on  the  east  side  of 
the  lake  and  have  focal  depths  of  less  than  8 km. 

Studies  of  Induced  earthquakes  from  reservoir  loading  require  a network 
of  seismographs,  adequate  for  the  approximate  location  of  small  local  earth- 
quakes, to  be  in  operation  before  the  Impounding  of  the  reservoir  (Bolt  and 
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Hudson,  1975).  Without  such  a network,  it  is  usually  impossible  to  establish 
the  seismicity  of  the  area  prior  to  closure.  Thus,  the  extent  to  which  local 
earthquakes  were  a consequence  of  the  reservoir  or  a part  of  the  more  general 
seismic  pattern  cannot  be  decided.  Such  a decision  is  essential  to  an  evalu- 
ation of  the  probable  location  and  size  of  future  shocks,  and  thus  is  of 
lassedi.\te  practical  Importance. 

1.3  Monitoring  of  Earthquake  Activity  near  Reservoirs 

In  the  preclosure  stage  of  a reservoir,  where  the  main  purpose  is  to 
establish  if  any  local  earthquakes  occur  usually  at  all,  a minimum  network 
of  three  to  five  short-period  vertical  component  seismometers  is  often 
recommended.  For  such  a network,  a rough  but  adequate  assessment  of  back- 
ground earthquake  frequency,  location  (using  P and  S waves),  and  magnitude 
can  be  made.  If  local  earthquakes  do  in  fact  occur,  the  network  can  be 
expanded  with  additional  seismometers  as  near  as  possible  to  the  active  area. 

After  closure,  the  usual  advice  is  to  operate  at  least  a four  to  six 
station  network  for  a period  extending  some  years  beyond  the  time  when 
maximum  impounding  is  complete.  If  a sequence  of  earthquakes  does  occur, 
then  the  network  is  usually  densified.  In  such  cases,  special  research  is 
often  warranted  and  careful  theoretical  considerations  must  be  brought  to 
bear  on  the  optimal  location  of  the  network.  Normally,  the  requirement  is 
to  obtain  an  accuracy  in  location  of  earthquake  foci  of  about  1 km  so  that 
correlations  with  geological  faults  can  be  made.  In  the  following  sections, 
an  account  is  given  on  how  such  a precision  can  be  obtained  with  local 
networks  of  various  sizes  and  configurations. 

In  practice  at  present,  sites  of  the  sensitive  seismographs  are  usually 
selected  based  on  practical  considerations  such  as  accessibility  and 
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avoidance  of  construction.  In  the  best  esses,  some  individual  site  selec- 
tion is  made  depending  on  cue  local  geological  structures.  The  instruments 
are  best  located  on  outcrops  of  basement  rock,  and  in  order  to  reduce  the 
microseismic  noise,  they  should  be  remote  as  can  be  achieved  from  construc- 
tion activities,  streams,  quarries,  spillways,  etc.  It  is  helpful  in  these 
respects  to  make  field  surveys  of  the  relative  background  microseismic  noise 
at  the  prospective  sites,  using  a portable  seismographic  recorder  before 
locations  are  finalized. 

There  are,  however,  more  fundamental  questions  concerning  the  optimal 
design  and  layout  of  the  network.  Let  us  consider  the  following  three  basic 
problems  to  be  solved  jointly:  a)  What  is  the  appropriate  velocity  struc- 
ture or  equivalently,  the  travel-time  curves  for  the  region  in  question? 

b)  How  can  the  specified  precision  of  hypocentral  location  be  obtained?  and 

c)  How  can  a minimer.  number  of  stations  be  involved?  Commonly,  at  present, 
these  questions  are  only  partially  treated.  First,  sometimes  chemical  ex- 
plosions are  fired,  perhaps  in  connection  with  quarrying  in  the  vicinity  of 
the  dam  construction,  to  obtain  observed  travel-time  curves  that  can  be  used 
in  the  earthquake  location  work.  More  usually,  already  available  travel-time 
curves  or  velocities  from  other  tectonic  regions  are  adopted  as  approximations 
by  analogy.  Secondly,  some  genetal  geometrical  considerations  are  used  to 
govern  the  network  configuration.  For  example,  sites  of  seismometers  are 
selected  to  be  as  uniformly  spread  in  azimuth  around  the  reservoir  as  possible. 
The  spacing  between  the  stations  is  kept  as  equal  as  practical,  and  is  usually 
selected  so  as  not  to  exceed  30  km  or  be  less  than  5 km. 

In  the  following  sections  the  three  Joint  problems  above  are  treated 
in  essentials.  In  the  course  of  this  discussion,  guidelines  •’■c  set  down 
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on  what  kind  of  seismic  readings  should  be  made  and  with  what  precision 
in  order  to  answer  questions  that  arise  concerning  seismicity  near  to  a 
reservoir.  As  a tool  in  the  analysis,  a simple  program  for  a high  speed 
computer  was  written  for  efficient  and  robust  numerical  estimation  of 
location  of  the  seismicity  near  to  a reservoir.  It  is  called  GHYP1. 


1 


2.0  THE  EARTHQUAKE  LOCATION  PROBLEM 


2.1  Crustal  Model  and  Travel  Times 

The  simplest  model  for  location  of  an  earthquake  focus  (l.e.,  hypocenter) 
Is  shown  In  Figure  1.  Four  parameters  need  to  be  determined  for  the  single 
focus  F,  namely  latitude,  longitude,  focal  depth,  and  origin  time.  The  epi- 
center Is  the  point  on  the  surface  directly  above  the  focus.  Stations  In 

* 

the  region  S^,  S2>  etc.  receive  P and  S waves  from  the  focus.  The  arrival 
times  of  these  waves  are  read  from  seismograms  and  used  to  develop  a system 
of  equations  for  the  determination  of  the  focal  parameters.  In  Figure  1, 
the  P and  S velocities  in  the  rocks,  called  and  V2,  are  constant,  so  that 
the  ray  paths  from  the  focus  to  the  stations  are  straight  lines.  For  a 


rectangular  set  of  spatial  axes,  the  coordinates  X,  Y,  and  Z of  a prelim- 
inary location  relative  to  a local  origin  are  related  by  Pythagoras'  theorem 


by  the  equation 

V(T  - Tq)  - (X2  + Y2  + Z2)**  (1) 

where  TQ  is  the  origin  time  and  T is  the  arrival  time  of  the  wave  in  question 
from  the  focus  to  the  station.  Call  this  case  Model  A. 

The  relation  between  the  coordinates  and  the  travel  time  is  thus  a non- 
linear one.  Rather  than  the  arrival  times  or  travel  times  of  the  waves,  the 
residual  times  (differences  between  the  observed  arrival  time  and  expected 
arrival  time  calculated  from  the  assumed  velocity  model)  are  used.  The  usual 
procedure  is  to  linearize  the  problem  on  the  assumption  that  the  adjustments 


For  an  elementary  treatment  of  seismic  waves  and  sources,  see  B.A.  Bolt, 
"Earthquakes  - A Primer,"  W.H.  Freeman,  1978. 
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AX,  AY,  AZ,  and  AT  to  some  provisional  location  X,  Y,  Z,  and  origin  time  T0 
are  small.  In  some  circumstances,  particularly  for  the  focal  depth  estimate 
Z,  this  is  not  the  case,  and  recently  methods  of  non-linear  regression  have 
been  used  (Marquardt,  1963). 

In  principle,  four  readings  of  arrival  times  of  P or  S waves  are 
sufficient  to  solve  for  the  four  unknowns,  AX,  AY,  AZ,  AT.  There  are, 
of  course,  measurement  errors  in  these  readings  from  seismograms  and  there 
are  variations  in  travel  time  due  to  changes  in  geological  structure  in 
various  directions.  Therefore,  more  readings  than  unknowns  are  obtained 
whenever  possible  and  the  resulting  overspecified  set  of  equations  solved 
by  least-squares.  A computer  program  GHYP1* based  on  this  simple  model  is 
described  in  the  following  sections. 

A more  complicated  Model  B,  meant  to  take  account  of  layering  in  the 
crustal  rocks,  is  shown  in  Figure  2 for  the  case  of  two  layers  over  a half 
space.  Here  waves  from  the  focus  can  reach  the  station  along  three  least 
time  paths.  The  P phases  in  this  case  are  denoted  as  Pj^,  P^  and  P^. 

The  first  path  is  direct,  but  the  second  two  are  diffracted  along  the 
boundary  between  the  layers.  Each  one  of  these  phases  could  arrive  first 
on  the  seismogram,  depending  upon  the  layer  thicknesses  and  the  distance 
between  the  station  and  hypocenter.  If  all  three  phases  can  be  read  on 
seismograms,  they  provide  additional  information  to  determine  the  four  focal 
parameters.  Widely  used  programs  of  this  type  of  which  the  model  is  composed 
of  discreet  layers  are  the  LOCAL  program  used  at  the  University  of  California, 
Berkeley  (Bolt  and  Turcotte,  1964)  and  a modificatioi  derived  from  it,  HYPO  71, 
used  at  the  U.S.  Geological  Survey  (Lee  and  Lahr,  IS. 3). 

'‘•This  program  is  available  from  the  Selamographic  Station,  University 

of  California,  Berkelay,  California  94720. 


A mathematical  advantage  of  this  type  of  model  can  be  seen  from  Figure 
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2.  The  upgoing  and  downgoing  rays  place  constraints  on  the  focal  depth 
from  above  and  below.  Thus,  if  there  are  stations  around  a reservoir  placed 
at  different  distances  from  the  focus,  the  first  arriving  and  waves 
(say)  will  give  a strong  estimate  of  focal  depth. 

The  layered  model  is  not  as  suited  for  study  of  localized  seismicity 
around  a reservoir  as  it  is  in  the  case  of  regional  earthquakes,  where 
stations  range  in  distance  from  tens  to  hundreds  of  kilometers  from  the  focus. 
In  the  case  of  reservoir  seismicity,  the  distances  between  the  focus  and  the 
station  are  usually  less  than  20  km,  so  that  no  seismic  phases  coming  from 
the  deeper  crust  and  being  diffracted  along  the  discontinuity  between  the 
crust  and  the  mantle  are  recorded  usually  on  the  seismograms.  It  is,  there- 
fore, unlikely  that  seismograms  will  show  discreet  wave  onsets  like  P^,  Pj^ 
and  P^  which  can  be  read. 

An  alternative,  called  Model  C,  is  to  introduce  a continuous  velocity 
increase  with  depth  in  the  upper  part  of  the  crust.  The  rays  in  this  case 
are  illustrated  in  Figure  3.  There  will  be  continuous  bending  of  the  rays 
for  laws  which  have  the  velocity  increasing  with  depth.  Mathematically,  the 
effect  on  the  solution  of  equation  (1)  is  similar  to  the  discreet  layer  case, 
because  the  rays  may  leave  the  focus  at  angles  either  above  or  below  the 
horizontal  through  that  point,  depending  on  the  distance  of  the  station  from 
the  epicenter.  A computer  program,  GHYP2,  based  on  this  model  has  been 
written  and  is  used  at  the  Seismographic  Stations,  University  of  California, 
Berkeley,  but  is  not  described  further  in  this  report. 
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2.2  Relative  Distribution  of  Hypocenters  and  Stations 


2.3  Errors  from  Timing,  Recording  and  Readlni 


Modern  networks  now  consist  of  Instrumentation  for  which  clock  errors 
ere  quite  smell.  Sometimes  s network  is  joined  by  s telemetry  to  a central 
recording  station,  and  in  this  case  relative  errors  are  essentially  zero. 
Even  if  each  station  has  its  own  crystal  clock,  however,  timing  errors  are 
not  likely  to  be  in  excess  of  0.1  second.  More  serious  are  the  errors  from 
the  variations  in  the  drum  rate  and  in  picking  the  onset  of  the  phase  in 
question.  Generally,  errors  up  to  i.3  sec  are  common,  and  when  the 
possibility  of  misinterpretation  of  phases  is  included,  larger  errors  in 
excess  of  half  a second  occur  with  a finite  frequency.  The  problem  of 
optimum  location  of  earthquakes  around  a reservoir  is  therefore  a statistical 
one,  and  the  analysis  of  the  problem  must  incorporate  the  statistical  vari- 
ations. 

Two  considerations  arise.  First,  the  purpose  of  making  a joint  location 
of  foci  of  several  earthquakes  rather  than  calculating  the  position  of  earth- 
quake foci  separately  (which  has  been  the  tradition  in  seismology  up  to  the 
present)  is  to  accumulate  a sample  of  readings  for  a particular  station.  In 
this  way  the  incorrect  readings  can  be  detected  against  the  mean  value  of  the 
sample.  Thus,  if,  say,  ten  local  earthquakes  are  located  together,  then  each 
station  of  the  network  will  have  approximately  ten  readings  of  the  P phase 
and  ten  readings  of  the  S phase.  These  will  provide  a stable  estimate  of  the 
mean  and  variance  of  the  station  measurements  arising  from  timing,  recording, 
and  reading  errors. 

Secondly,  the  measurement  errors  must  be  treated  in  the  program  by  means 
of  a suitable  weighting  scheme.  The  computational  algorithm  systematically 
tries  to  improve  the  location  estimates  by  successive  iterative  adjustments. 
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A necessary  condition  for  solving  such  an  inverse  problem  is  that  the 
forward  problem  is  tractable.  As  we  pointed  out  in  section  2.0,  the  forward 
problem  is,  in  fact,  a non-linear  problem.  It  can  be  treated  as  such  or 
made  more  tractable  by  linearizing  it.  The  latter  is  the  approach  we  take 
here. 

As  argued  in  the  last  section,  it  is  statistically  advantageous  not 
to  locate  earthquakes  one  by  one,  but  rather  accumulate  the  times  from  a 
number  of  earthquakes  near  the  reservoir  and  then  locate  them  all  simultane- 
ously in  the  one  joint  least-square  computation  (Douglas,  1967;  Dewey,  1971; 
Crosson,  1976).  Let  us  suppose  that  there  are  n earthquakes  and  m sta- 
tions in  the  network.  Suppose  further  that  the  number  of  model  parameters 
to  be  adjusted  (seismic  velocities)  is  it.  The  problem  then  is  one  of  deter- 
mining up  to  An-fi-Hn  unknowns.  There  are  n x m observations,  if  ona  phase 
is  read  at  the  station,  or  n x 2m  observations,  if  both  P and  S arc  read  at 
one  of  the  stations.  In  the  latter  case,  for  example,  there  may  be  five 
earthquakes  recorded  at  six  stations  and  a P and  an  S velocity  correction 
to  be  estimated.  Thus  there  would  be  28  unknowns  to  be  determined  from  60 
observations.  Such  a problem  is  soluble  in  the  least  squares  sense. 

Prom  an  analytic  point  of  view,  the  problem  can  be  thought  of  as  one 
in  the  analysis  of  variance.  This  statistical  method  is  well-worked  out 
and  provides  a valuable  theoretical  framework  for  analysis. 
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3.0  MATHEMATICAL  FORMULATION 

3.1  The  Algorithm  for  Simultaneous  Inversion 

We  now  set  up  the  general  problem  in  mathematical  form  along  the 
llnea  explained  In  2.0  for  Models  A and  B. 

Let  the  suffix  1 refer  to  the  ith  earthquake  (1-1,  ...  , n)  and  j 
refer  to  the  jth  station  (J  - 1,  ...  , m) . Let  the  time  adjustment  special 
to  the  Jth  station  be  a ^ . Consider  the  1th  earthquake  being  recorded  at  the 
jth  station.  For  simplicity,  wo  consider  only  the  first  P wave  to  be 
measured.  Equations  for  the  additional  phases  ( S , for  example)  can  be 
easily  added.  Then  the  residual  (observed  minus  calculated)  In  the  measured 
P arrival  time  is  r^  and  this  Is  associated  with  a random  error  term  . 
Then  we  may  write,  correct  to  the  first  order, 

rlj  + clj  * Atl  + LAX1  + **^1  + NAZi  + HAV1  + A«j  » (2) 

where 
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These  partial  derivatives  must  be  calculated  from  the  appropriate  formula 
linking  travel  time  with  velocity  and  distance  along  a ray.  In  the  case  of 
the  constant  velocity  model,  the  formula  is  equation  (1). 

There  will  be  an  equation  like  (2)  for  each  observed  phase  at  each 
station  for  each  earthquake.  This  set  of  linear  equations  of  condition  can 
then  be  solved  for  the  unknowns  Atj,  AXj,  AY ^ , AZ^,  AV^,  and  Aa^  (all  1 and 
J).  The  solution  Involves  the  Inversion  of  a matrix  (see  below)  and  often 
has  numerical  problems  associated  with  ill  conditioning.  The  aim  of  network 
design  is  to  remove  this  111  conditioning  as  much  as  possible. 
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The  set  of  linear  equations  of  condition  may  be  written  as 

Ax  - b,  (3) 

where  x are  the  unknown  adjustments  and  b are  the  residuals. 

Each  equation  can  be  weighted,  by  use  of  a weighting  matrix  W,  to 

allow  for  relative  uncertainties  in  the  observations.  In  general,  there 

T 

are  more  observations  than  unknowns  so  that  A A is  nonsingular.  Therefore, 

T -1  T 

x - (A  A)  A A b . (4) 

When  there  are  fewer  observations  than  unknowns,  the  generalized  least- 
squares  inverse  may  be  used  (Bolt,  1970), 

x - AT(AAT)-1b  . (5) 

In  addition  to  the  condition  equations,  a number  of  linear  constraints 
must  be  imposed  in  general  (such  as  requiring  the  sum  of  station  adjustments 
to  be  zero).  These  define  the  constraint  matrix 

Cx  - d . (6) 

We  then  solve  in  the  computer 
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for  the  hypocentral,  velocity,  and  station  adjustments  x and  the  Lagrange 
multipliers  X (not  used) . 

In  addition,  the  allowable  range  of  adjustments  of  the  station  terms 
Aaj  and  P velocity  AV^  was  controlled  by  augmenting  the  matrix  equation 
(2)  with  penalty  equations  Aa^  ■ 0,  AV^  ■ 0.  Deviations  from  these  values 
are  penalized  but  not  prevented  by  applying  appropriate  weights. 
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It  will  be  clear  from  (3)  and  (2)  Chat  the  matrix  A will  lead  to  111- 

T 

conditioned  lnveraes  (l.e.,  a nearly  singular  A A)  If  the  equations  of 
condition  are  almost  parallel.  This  problem  may  be  reduced  by  ensuring 
that  the  ratios  of  the  coefficients  L,  M,  N,  H of  (2)  are  as  different  as 
possible  for  each  earthquake-station  pair.  That  Is,  the  stations  should  not 
be  at  the  same  distance  from  the  reservoir.  In  addition,  the  coefficients 
for  S waves  usually  differ  considerably  from  those  for  P waves,  so  that 
both  P and  S wave  onsets  should  be  read  wherever  possible. 

The  solution  Is  programmed  to  proceed  by  a series  of  steps.  At  each 
stage,  the  corrections  are  reapplied  to  the  previous  set  of  values  of  the 
hypocentral  coordinates,  origin-times,  velocities,  and  station  adjustments. 

This  yields  a new  set  of  residuals  b and  a new  set  of  coefficients  A. 

Iterations  proceed  until  the  values  of  x become  less  than  a pre-specified 
limit. 

In  the  case  of  GHYP1,  which  uses  (4),  the  weight  matrix  W is  automatically 
calculated  at  each  iteration  from  a Pearson  Type  VII  law  whose  w dth  is  pro- 
portional to  the  variance  of  the  residuals  at  each  stage.  Thus,  as  the 
solution  converges,  the  weighting  becomes  more  severe. 


3.2  Optimization  of  Station  Distribution 

Theoretically,  the  properties  of  the  matrices  of  condition  may  be 
examined  to  determine  their  resolving  power.  Usually  there  is  partial 
redundancy  in  the  constraints  imposed  by  the  measurements  available.  For 
example,  two  stations  close  together  will  provide  two  equations  which  are 
almost  parallel.  This  ill-conditioning  shows  up  by  examining  the  relative 
magnitudes  of  the  eigenvalues  of  the  matrices  involved. 

For  the  present  introductory  summary,  a numerical  trial-and-error 
scheme  has  been  used  to  illustrate  the  effect  of  varying  the  station  lo- 
cations around  a small  region  such  as  a reservoir. 

Consider  16  stations  distributed  around  the  reservoir,  as  in  Figure  A, 
in  two  concentric  circles , radii  5 and  10  km.  Stations  and  Slg  are  in 
the  inner  ring  and  stations  to  S£g  are  in  the  outer  ring.  The  epicenter 
of  the  earthquake  is  assumed  to  be  at  the  center  of  the  circle. 

The  aim  of  the  experiment  was  to  investigate  the  effect  of  simple  changes 
in  the  arrangements  of  these  16  stations.  Crustal  Model  B was  used  with  a 1 
km  thick  surficial  layer  (V^  » 2.0  km/sec)  and  a 2 km  thick  intermediate  layer 
(Vj  * 3.0  km/sec).  The  P velocity  in  the  under  layer  was  “ 5.5  km/sec. 

For  a focal  depth  of  2 km,  the  travel  time  (errorless)  to  the  inner  ring 
stations  was  2.08  sec  along  a direct  path  (Pjj)  an<*  t*ie  errorless  time  to  the 
outer  ring  stations  was  3.13  sec  along  the  refracted  path  (P23  ~ see  Figure  2). 
By  the  throw  of  a die  these  times  were  then  varied  for  each  station  by  a 
random  station  adjustment  ±0.3,  ±0.2,  ±0.1.  The  resulting  adjusted  travel 
times  were  then  used  with  program  LOCAL  to  locate  the  epicenter  of  the 
earthquake.  (Focal  depth  was  held  fixed.) 
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Nine  cases  were  run  with  various  selections  of  stations.  For  example, 
Case  1 had  all  8 stations  in  the  inner  ring  and  none  in  the  outer  ring; 

Case  9 had  S^,  S^,  s^g»  S27»  an<*  S28  olnitte^«  i.e.,  all  stations  in  the 
northeast  quadrant  absent. 

Details  of  the  results  will  not  be  given  here,  but  two  results  illus- 
trate the  analysis. 

(a)  There  is  a dramatic  loss  of  resolution  when  all  stations  in  one 
quadrant  were  omitted  (Case  9).  Compared  with  the  solution  with 
all  16  stations  present,  the  standard  error  of  one  observation 
changed  from  0.31  sec  to  0.37  sec  and  the  latitude  and  longitude 
by  0.4  km. 

(b)  Inclusion  of  the  second  ring  does  not  significantly  alter  the 
resolution  of  the  epicenter  when  the  full  equispaced  inner  ring 
is  in  operation  but  drastically  improves  the  estimation  of  focal 
depth. 

We  now  illustrate  the  ideas  developed  in  the  previous  sections  by 
discussing  two  examples  of  sequences  of  earthquakes  near  reservoirs  that 
occurred  when  only  local  networks  were  available  for  location.  One  is  a 
sequence  near  the  Oroville  reservoir  in  central  California,  and  the  other 
a swarm  of  small  earthquakes  near  the  Briones  reservoir  in  the  hills 
eaet  of  San  Francisco  Bay.  Neither  network  had  the  optimal  geometrical 
arrangement  shown  in  Figure  4,  but  where  possible  this  should  be  approximated 
as  closely  as  practical. 
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4.0  CASE  STUDIES 

4.1  The  Orovllle  Sequence.  August  1975 

The  main  shock  (magnitude  5.7)  of  the  sequence  (see  Figure  5)  occurred 
in  the  afternoon  of  August  1,  1975,  at  1320  PDT.  The  region  Is  one  of  low 
seismicity  with  generally  a few  minor  earthquakes  a year  within  50  km  of 
the  dam  site.  Nevertheless,  before  dam  construction,  seismographs  were 
installed  in  1963  about  a kilometer  north  of  the  dam  in  order  to  monitor 

1 i 

the  background  seismicity.  These  instruments  detected  no  change  in  the  low 
level  of  earthquake  occurrence  within  30  km  of  the  reservoir  during  filling 
or  after  the  reservoir  was  raised  to  its  highest  elevation  in  1969. 

On  June  28,  1975,  a few  small  shocks  were  recorded  southwest  of  the 
Oroville  reservoir.  Some  additional  portable  seismographic  stations  were 
installed  to  keep  better  track  of  the  position  of  the  earthquakes.  About  20 
small  shocks  were  recorded  through  July  in  the  same  general  area,  the  largest 
of  which  was  magnitude  4.7.  After  the  main  shock  on  August  1,  many  after- 
shocks followed  in  the  subsequent  weeks  and  many  more  temporary  seismographs 
were  operated  in  a dense  local  network.  Locations  of  51  earthquakes  in  the 
sequence  are  plotted  in  Figure  5.  The  main  shock  was  located  about  10  km  to 
the  south  of  Lake  Oroville  and  the  foreshocks  and  aftershocks  define  a region 
of  area  10  km  by  14  km  to  the  south  of  the  dam  (Morrison  et  al.  1976).  The 
focal  depths  ranged  from  10  km  to  the  west  of  this  zone  to  near  surface 
values  to  the  east. 

The  hypocenters  in  Figure  5 were  calculated  one  by  one,  using  the  routine 
program  LOCAL,  using  P readings  from  nine  stations  with  epicentral  distances 
from  7 to  195  km.  Only  5 stations  were  available  at  the  time  of  the  main 
shock  at  distances  of  40  km  or  less  from  the  epicenter.  Thus,  these  initial 
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hypocenter  estimates  do  not  come  from  the  most  satisfactory  station  ar- 
rangement. Nevertheless,  a sub-set  of  the  aftershocks,  when  temporary 
stations  were  available,  provides  a useful  data  set  to  test  and  demonstrate 
the  use  of  the  Joint  location  program  GHYP1  (see  Appendix). 

In  this  test,  9 near  stations  (distances  from  2 to  25  km)  were  selected 
with  even  acimuthal  distribution  about  9 aftershock  hypocenters.  There  were 
P and  S readings  at  each  station  for  each  shock  yielding  162  observations 
and  47  parameters  to  be  estimated  (see  section  3).  Initial  values  for  V^, 

P,  and  S velocities  were  ■ 5.9  km/sec  and  - 3.47  ka/sec,  respectively. 

After  10  iterations  with  GHYP1,  the  final  solution  Is  summarised  in 
Table  1.  The  table  also  shows  the  solutions  at  the  first,  third,  and  fifth 
iterations.  (The  standard  errors  of  the  estimates  are  given  In  parentheses.) 
It  Is  of  interest  to  compare  the  final  solution  using  the  Joint  location 
method  with  the  earlier  independent  solutiona  (last  column  GS)  obtained 
with  a wider  station  coverage  but  an  individual  location  method.  Epicenters 
have  moved  by  a kilometer  or  so  relative  to  each  other  while  focal  depths 
have  changed  by  up  to  3 km.  The  station  time  adjustments  are  small  for 
both  P and  S,  the  largest  being  about  0.1  sec.  The  solution  suggests 
Vj  ■ 6.07  i 0.06  km/sec  and  Vj  ■ 3.49  i 0.02  km/sec  are  Improved  values  for 
the  mean  velocities. 

As  an  illustration,  the  computer  output  for  iteration  4 and  the  magnitude 
1 aftershock  on  14  July  1976  is  shown  in  Figure  6 (see  Appendix).  The  first 
column  of  residuals  are  from  P waves,  the  second  column  are  from  S waves.  The 
upper  triangle  of  the  correlation  matrix  is  tabulated  in  each  case  to  indicate 


the  amount  of  correlation  between  the  estimated  adjustments  to  latitude, 
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TABU  1 


GROUP  LOCATION  OF  OROVILLE  EARTHQUAKES 
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OT 
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18.041.026) 
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17.84 

♦ 

28.14 

29 . 76(  , 11) 
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29.561.05) 

29.65 

X 
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29.73 

h 
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4 . 54( , 34) 

4.631.17) 
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OT 

26.90 
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27 .521.040) 
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28.14 

23 . 80( ,16) 

23.861.09) 

23.891.08) 

24.19 

X 

29.60 

29.06( , 22) 

28.961.13) 

28.981.11) 

28.98 

h 

6.73 

6.51(.46) 

6.471.24) 

5.831.25) 

5.66 

1003Febl9 

OT 

09.45 

09. 59( .050) 

09.591.025) 

09.691.035) 

09.48 

♦ 

28.14 
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30.221.08) 

30.161.07) 

30.31 

X 

29.60 

31 . 16( . 26) 

31.421.15) 

31.481.13) 

31.43 

h 

6.73 

8. 53( . 28) 

8.691.14) 

8.351.15) 

7.91 

1214F*bl9 

OT 

18.38 

18.65( .041) 

18.651.021) 
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♦ 

28.14 

24.94( . 12) 

24.951.07) 

24.961.07) 

X 

29.60 

29.341.20) 

29.341.12) 

29.391.10) 

h 

6.73 

7.64(.25) 
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OT 
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23.961.040) 
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X 
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30.001.16) 

30.081.13) 
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h 

6.73 

4.861.46) 
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3.751.27) 

3.65 
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OT 

34.40 
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35.12 
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28.14 
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27.191.07) 

27.161.06) 

27.23 

X 

29.60 

30.081.20) 

30.251.11) 

30.431.09) 

30.36 
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6.73 

4.051.36) 

3.781.19) 

3.431.18) 

3.57 

0007Marl3 

OT 

44.53 

44.471.053) 

44.491.027) 

44.601.038) 

44.37 

4 

28.14 

26.351.14) 

26.381.08) 

26.381.07) 

26.62 

X 

29.60 

32.371.22) 

32.301.13) 

32.251.12) 

31.93 

h 

6.73 

9.141.29) 

9.061.15) 

8.681.16) 

8.78 

030lJull4 

OT 

57.04 

57.431.041) 

57.441.021) 

57.551.033) 

57.30 

4 

28.14 

24.311.11) 

24.321.07) 

24.341.07) 

24.55 

X 

29.60 

30.751.17) 

30.781.11) 

30.801.09) 

30.70 

h 

6.73 

6.191.25) 

6.101.13) 

5.581.16) 

5.90 

0304Jull4 

OT 

14.47 

14.771.050) 

14.841.025) 

14.941.033) 

14.72 

4 

28.14 

24.301.12) 

24.331.07) 

24.341.07) 

24.56 

X 

29.60 

30.431.22) 

30.721.12) 

30.751.10) 

10.68 

h 

6.73 

6.701.29) 

6.311.15) 

5.811.16) 

5.97 

ICont (mini) 
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TABLE  1 (Concluded) 


INITIAL 

FINAL 

PARAMETER 

SOLUTION 

3RD  ITERATION 

5TH  ITERATION 

SOLUTION 

CS»* 

P ADJUSTMENT 

TAB 

0 

.025(.011) 

.023(.010) 

.03 

COX 

0 

-.047( .013) 

-.048( . Oil) 

.06 

LUV 

0 

,042(.010) 

•048( .009) 

.04 

ORV 

0 

,033(.013) 

. 043( .01 1) 

-.03 

KPK 

0 

-,048( .014) 

-.064( .015) 

.07 

FIG 

0 

-,013( .010) 

-.014 (.010) 

WYN 

0 

,012( .009) 

.014(.010) 

.05 

LON 

0 

-,055(.011) 

-.059( .010) 

.08 

HON 

0 

.051(.013) 

. 058( . Oil) 

-.08 

S ADJUSTMENT 

TAB 

0 

•031( .014) 

.049( .013) 

-.05 

COX 

0 

-.022( .016) 

-.044(.014) 

-.12 

LUV 

0 

.016( .012) 

.025(.011) 

-.14 

ORV 

0 

. 140( .014) 

. 161 ( .012) 

-.08 

KPK 

0 

-.114( .019) 

-.093( .017) 

-.19 

FIG 

0 

-.001(.014) 

-.023( .014) 

WYN 

0 

,002( .012) 

-.010( .013) 

.05 

LON 

0 

-.126(.015) 

-. 149( .013) 

HON 

0 

,072( .016) 

,085( .015) 

VELOCITIES 

P 

5.9 

6.07( .057) 

5.9 

S 

3.47 

3.49( .016) 

3.4 

USGS  unpublished  data. 


longitude,  depth,  etc.  The  relatively  large  1-4  element  indicates  a 
high  correlation  between  adjustments  in  origin  time  and  focal  depth,  while 
the  smaller  2-3  element  indicates  low  correlation  between  adjustments  in 
latitude  and  longitude.  This  matrix  can  be  used  to  calculate  uncertainty 
ellipses  if  required  (see  Dewey,  1971).  The  azimuthal  distribution  (epicenter 
to  station  ) can  be  used  to  compare  with  Figure  4. 

4.2  The  Briones  Hills  Swarm,  January  1977 

A swarm  of  small  to  moderate  shallow-focus  earthquakes  (Bolt  et  al, 

1977)  occurred  in  the  vicinity  of  Briones  Hills  in  Contra  Costa  County  on 
the  weekend  of  January  7 to  10,  1977  (local  time).  In  four  days,  69  earth- 
quakes of  magnitude  1.0  or  greater  had  been  recorded  at  the  University  of 
California  station  BKS  situated  10  km  away.  The  largest  earthquake  had 
magnitude  4.3.  The  epicenters  were  located  near  the  Briones  and  San  Pablo 
reservoirs  and  there  was  considerable  interest  in  a possible  link  between 
the  swarm  and  reservoir  loading  (see  Figure  7). 

High-quality  hypocentral  parameters  could  be  estimated  for  eight  of 
the  earthquakes  using  P readings  and  six  nearby  stations  (4  km  < distance 
< 23  km)  and  S readings  from  the  BKS  recordings.  The  epicenters,  calcu- 

lated separately  for  each  earthquake  using  Model  B and  the  program  LOCAL, 
are  plotted  in  Figure  7 as  crosses.  Bars  denote  standard  errors  in  the  co- 
ordinates of  the  epicenters.  Fecal  depths  ranged  from  7 km  to  12  km.  The 
best  locations  indicated  that  the  earthquake  sources  were  not  under  Briones 
reservoir  and  probably  not  associated  with  it. 

This  set  of  earthquakes  provides  a further  suitable  test  of  the  joint 
location  scheme  discussed  earlier.  The  P and  S times  for  the  eight  earthquakes 
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TABLE  2 

GROUP  LOCATION  OF  BRIONES  HILLS  EARTHQUAKES 


♦ - 37*  + N X - 122*  + W 


PARAMETER 

06SSJan08 

OT 

INITIAL 

SOLUTION 

3RD  ITERATION 

5TH  ITERATION 

FINAL 

SOLUTION 

PAPER* 

50.84 

50.23(.301) 

50.20( .126) 

50.48( .170) 

50. 55( • 28) 

♦ 

54.86 

54.04(1.18) 

53.98(.53) 

53. 64( .41) 

54 ,19( >94) 

X 

9.17 

11.49(1.11) 

11. 55( . 50) 

12.1K.44) 

11.09( .88) 

h 

7.64 

13.98(1.82) 

14.02(.77) 

13.04(.96) 

11.74(1.82) 

0717Jan08 

OT 

33.79 

33.72(.278) 

33.67(. 116) 

33.98(.136) 

33.94(.09) 

♦ 

54.86 

53.70(1.08) 

53.62(.48) 

53.45(.37) 

53.80(.32) 

* 

9.17 

11.57(1.02) 

11.69( .46) 

11.93( .37) 

11.31(.30) 

h 

7.64 

10.98(1.76) 

11.21(.75) 

10.03( .81) 

9. 33( . 64) 

0858Jan08 

OT 

13.79 

13.63(.282) 

13.57( .118) 

13.90( .142) 

13.87( .10) 

♦ 

54.86 

53.77(1.08) 

53.69( .49) 

53. 53( .37) 

53.88(.35) 

X 

9.17 

11.32(1.01) 

11.43( .46) 

11.66(.32) 

11.05(.32) 

h 

7.64 

11.44(1.79) 

11 . 69( . 76) 

10.46( .85) 

9.65( .79) 

0938Jan08 

OT 

07.49 

07.30(.277) 

07. 22( . 116) 

07.53(.141) 

07.49( .11) 

♦ 

54.86 

54.26(1.04) 

54 . 16( .47) 

54.00( .35) 

54.13(.37) 

X 

9.17 

11.14( .97) 

11.3K.44) 

11.56(.36) 

10.97( .34) 

h 

7.64 

10.94(1.84) 

11.38(.78) 

10. 18( .86) 

9.47( . 79) 

0941Jan08 

OT 

02.69 

02.50( .283) 

02.43( .120) 

02. 64( . 148) 

02.69( .15) 

♦ 

54.86 

54.27(1.07) 

54.06(.49) 

53.59(.40) 

54 .74(  .49) 

X 

9.17 

12.07(1.08) 

12.30( .50) 

13.24( .48) 

11 . 35( .48) 

h 

7.64 

11.25(1.82) 

11.61(.76) 

11.20( .81) 

9.34(1.09) 

0951Jan08 

OT 

55.19 

55.43(.259) 

55.32(.110) 

55.61(.125) 

55.57( .13) 

♦ 

54.86 

54. 54(.94) 

54.39(.43) 

54.27(.32) 

54 . 50( .45) 

X 

9.17 

11.36(.91) 

11.58(.42) 

11.79(.34) 

11.25( .43) 

h 

7.64 

8.80(1.87) 

9.51(.79) 

8.39( .81) 

7.61(1.07) 

0534Jan09 

OT 

16.69 

16.54(.226) 

16.49(.115) 

16.80( .135) 

16.71(.08) 

♦ 

54.86 

53.24(1.10) 

53.15(.49) 

52. 97 (.38) 

53.31(.29) 

X 

9.17 

11.18( .98) 

11.3K.44) 

11.55(.35) 

11 .06( , 26) 

h 

7.64 

10.89(1.75) 

11 . 12( .75) 

9.93(.81) 

9.44( .56) 

054  6 Jan  09 

OT 

40.19 

40.30(.262) 

40. 22( .110) 

40. 50( .124) 

40. 44 ( ,08) 

4 

54.86 

53.06(1.05) 

52.94(.47) 

52. 76( .37) 

53.15(.28) 

* 

9.17 

10.80( .89) 

10.99( .41) 

11 . 22 ( . 32) 

10.70(.24) 

h 

7.64 

9.38(1.77) 

9.79( .45) 

8.68(.78) 

8.12(.57) 

P ADJUSTMENT 

BKS 

0 

-.011(.026) 

.001 (.021) 

-.03 

BWR 

0 

-.257(.035) 

-.293(.032) 

-.19 

BRK 

0 

. 112( .034) 

. 170( .024) 

.06 

DSR 

0 

-.130( .038) 

-.135( .027) 

-.15 

BOL 

0 

-.006 (.043) 

-.010( .032 

-.04 

DUR 

0 

.291(.035) 

.267( .031) 

.25 

(Continued) 
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TABLE  2 (Concluded) 


PARAMETER 


INITIAL 

SOLUTION  3RD  ITERATION  5TH  ITERATION 


FINAL 

SOLUTION  PAPER* 


S ADJUSTMENT 
BKS 


VELOCITY 

P 

S 


5.6 

3.3 


Bolt  at  . 1977. 


. 127( .053) 


. 1 2 5 ( .039)  .23 


5.98( .144) 
3.23( .092) 


5.6 

3.3 
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5.0  CONCLUSIONS 


The  following  set  of  guidelines  for  the  design  of  a limited  network 


of  seismographic  stations  around  a reservoir  or  similar  facility  resulted 


from  the  studies  outlined  above 


1.  Location  of  stations 


a.  Survey  prospective  sites  for  microseismic  noise  level;  sites  with 


the  lowest  noise  levels  are  preferable 


fe.  Outcrops  of  basement  rock  are  preferable  sites  to  those  on  loose 


alluvium  or  fill 


c.  Spacing  between  adjacent  stations  should  be  approximately  equal 


and  should  not  exceed  30  km  or  be  less  than  5 km 


d.  The  station  network  should  be  approximately  centered  on  the  facility 


e.  Stations  should  be  located  around  the  facility  at  approximately 


equal  Intervals  of  azimuth 


2.  Number  of  stations 


a.  Four  stations  is  the  minimum  number  required  for  reasonably  precise 


location  if  both  S and  P phases  can  be  read  at  all  stations  for 


most  of  the  smaller  earthquakes 


b.  Seven  stations  is  the  recommended  minimum  number  if  an  S phase  can 


be  read  at  only  one  of  the  stations 


3.  Recommended  distribution  of  minimum  number  of  stations  (assuming 


region  of  interest  is  approximately  10  km  across) 


a.  Four  stations:  equilateral  triangle  15  km  on  a side,  with  the 


fourth  station  near  the  center 


b.  Seven  stations:  hexagon  10  km  on  a side  with  the  seventh  station 


near  the  center 


4.  Ocher  criteria 


a*  Onset  times  of  P and  S waves  should  be  measured  to  an  accuracy 
of  at  least  0. OS  seconds  and  preferably  to  0.02  seconds. 

b.  Station  locations  should  be  known  to  within  10  meters. 

c.  Magnification  of  each  station  should  be  set  to  a level  where  the 
excursions  of  the  trace  on  the  recording  device  due  to  the  back- 
ground microseisms  are  about  0.5  to  It  of  full  scale. 

When  P and  S measurements  come  from  a network  designed  using  the  above 
guidelines,  the  group  hypocenter  location  program  GHYP1  results  in  formal 
estimates  of  standard  errors  of  about  0.2  km  for  epicenter  location  and 
0.3  km  for  depth  (or  better).  These  uncertainties  are  sufficiently  small 
to  permit  correlation  with  surface  geological  features. 
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APPENDIX  A:  GHYP1  USERS  GUIDE 


Program  GHYP1  is  designed  to  locate  local  earthquake  sources  by 
groups,  rather  than  individually,  using  a small  array  (less  than  40  km, 
say)  of  sei8mographic  stations  and  a half-space  velocity  model  (see 
Figure  1) . GHYP1  can  accept  data  for  up  to  10  earthquakes  recorded  by 
up  to  10  stations  and  each  station  can  have  a P and/or  an  S observation 
for  each  earthquake. 

The  program  simultaneously  estimates: 

1.  hypocentral  parameters, 

2.  station  adjustments  (to  be  added  to  calculated  times),  and 

3.  P and  S propagation  velocities. 

The  estimation  is  made  by  first-order  adjustments  to  a starting 
solution  that  is  automatically  provided;  velocities  must  be  specified. 

A FORTRAN  listing  and  card  deck  of  the  program  GHYP1  are  available  at 
the  U.S.  Army  Engineer  Waterways  Experiment  Station  and  at  the  Seismo- 
graphlc  Station,  University  of  California,  Berkeley. 
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1.  Read  VP,VS,SMAX,NIT1,NIT2,NITMAX, (JS(I) ,1-1,5) 


Format  (3F10.2,8I5) 
where : 

VP  ■ initial  P wave  velocity  (km/ sec)  (5.6,  say) 

VS  ■ Initial  S wave  velocity  (km/sec)  (3.3,  say) 

SMAX  * maximum  step  size  (used  to  determine  if  solution  for  system 
of  equations  has  converged) . Suggest  SMAX  - VP*AT  where  AT 
is  resolution  in  reading  onset  times.  (SMAX  - 5.6x.02  - .1, 
say) 

NIT1  “ Number  of  iterations  solving  for  hypocentral  parameters  only, 
(usually  ■ 2 or  3) 

NIT2  “ Number  of  iterations  solving  for  hypocentral  parameters  and 
station  adjustments,  (usually  ■ 1,  2,  or  3) 

NITMAX  ~ Maximum  number  of  iterations.  Number  of  iterations  solving 

for  hypocentral  parameters,  station  adjustments,  and  velocity 
model  - NITMAX-NIT1-NIT2 . (NITMAX  - 7 to  12,  say) 

JS(1)  to  JS(5)  are  printing  sense  switches. 

JS(1)“0:  print  station  location  and  phase  data. 

JS(2)“1:  print  station  coordinates  and  parameters  used  in  setting 
up  the  size  of  the  equation  of  condition  matrix. 

JS(3)“1:  print  equations  of  condition. 

JS(4)“1:  print  the  variables  involved  in  solving  the  system  of  normal 
equations . 

JS(5)"1:  print  perturbation  of  parameters 
Note:  JS(2)  to  JS(5)  apply  to  each  iteration. 


2.  Read  In  station  data  (up  to  10  stations) 

Read  SNAME(I) ,SLATD(I) ,SLATM(I) ,SL0NGD(l) ,SL0NGM(I) 

Format  (A3,1X,F2.0,F5.2,1X,F3.0,F5.2) 
where : 

SNAME(I)  - station  code  (up  to  3 characters) 

SLATD(I)  - station  latitude  (degrees) 

I 1 

SLATM(I)  - station  latitude  (minutes) 

SLONGD(I)  ■ station  longitude  (degrees) 

SLONGM(I)  - station  longitude  (minutes) 

Note:  latitude  is  assumed  to  be  North  and  longitude  to  be  West. 

Repeat  for  each  station  and  after  the  last  station  place  and  end-of-record 
(BOR)  in  the  card  deck. 

3.  Read  in  phase  onset  time  data  (up  to  10  earthquakes). 

a.  Read  (ID(I) ,1-1,8)  - arbitrary  title  to  identify  the  event. 

Format  (8a10) 

b.  Read  STN(I) ,PH(I) ,PM(I) ,PS(I) ,SH(1) ,SM(I) ,SS(I) 

Format  (A3,1X,2F2.0,F5.2,1X,2F2.0,F5.2) 
where : 

STN(I)  - station  code 

PH, PM, PS  - P wave  onset  time  in  hours,  minutes,  and  seconds 
(to  .01  sec). 

SH,SM,SS  - S wave  onset  time  in  hours,  minutes,  and  seconds 
(to  .01  sec) 

Repeat  step  b for  each  station  observations. 

c.  BOR 

Repeat  steps  a,  b,  and  c for  each  event. 

A.  BOR  to  indicate  end  of  input  data. 
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Computational  Method 


GHYP1  uses  constrained  linear  least-squares  with  penalty  functions 
to  solve  a system  of  equations  and  the  techniques  of  analysis  of  variance 
are  used  to  estimate  the  standard  errors  of,  and  correlations  between, 
the  unknown  parameters  from  the  standard  error  of  the  residuals.  The 
observations  are  weighted  by  their  residuals  using  a Pearson's  Type  VII 
distribution. 

Constraint  equations  are  used  to  constrain  the  sum  of  the  P station 
adjustments  and,  if  applicable,  the  S station  adjustments, to  be  zero. 
These  constraints  are  necessary  to  avoid  a singular  system  of  equations 
due  to  unity  correlation  between  the  station  adjustments  and  the  origin 
times  of  the  earthquakes. 

Penalty  functions  are  used  to  restrain  the  perturbations  in  the 
station  adjustments,  the  velocities,  and  the  velocity  ratio  to  be  small 
(less  than  0.1,  say)  in  order  that  the  solution  to  the  system  of  normal 
equations  remains  stable. 

The  initial  location  for  each  epicenter  is  the  geometrical  center  of 
the  array  and  the  depth  is  0.7  of  the  average  radius  to  the  stations. 

The  initial  origin  time  depends  upon  the  initial  velocity  model. 


The  program  prints  out,  for  each  iteration  where  applicable: 
1.  The  hypocenter  for  each  earthquake,  including 


g.  standard  errors  of  parameters 


b.  correlation  matrix  (used  to  construct  error  ellipsoid  (uses 
left-hand  coordinate  system)) 

c.  perturbation  of  each  parameter  from  previous  iteration. 

2.  Station  data  for  each  earthquake,  including: 

§.  delta  (km) 

]}.  azimuth  (epicenter  to  station) 
c.  P onset  time,  weight,  and  residual 
£.  S onset  time,  weight,  and  residual. 

3.  Station  adjustments  and  standard  errors  for  P and  S. 

4.  P and/or  S velocities  and  standard  errors. 

Included  in  the  output  for  each  Iteration  is  an  estimate  of  the  number 


of  significant  figures  in  the  solution,  and  the  standard  error  of  the 
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In  accordance  with  letter  from  DAEN-RDC,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 
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